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ABSTRACT 


STXCAP (Stiff Circuit Analysis Program) is a FORTRAN IV, Version 
2.3, con 5 >uter program written for the CDC-6400-6600 computer series 
and SCOPE 3.0 operating system. It provides the circuit analyst a 
tool for automatically computing tlie transient responses and fre** 
quency responses of large linear time Invariant net^^orks,. both stiff 
and non~stiff. The circuit description and user's program input 
language is engineer-oriented, making single the task of using the 
program. 

Three volumes of documentation are available for the STIC/IP 
program; a theory manual, a user's manual, and a system's programmers 
manual. Volume I describes tlie engineering theories underlying 
STICAP and gives further references to the literature. Volume II, 
the user's manual, explains user interaction vjith the program and 
gives results of typical circuit design applications. Volume III 
depicts the program structure from a system's programmers vier^point 
emd contadns flow charts and other software doctmentation. 
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CHAPTER I 

* . 

IHTRODUCTIO:! 

Most computer aided netv;oric analysis programs are designed 
to relieve the user of the follo\7ing burdens: 

(a) Circuit translation - the «.btaining of tlie differential 
and algebraic equations governing tiie netv/ork, starting 
from an easily specified description of the circuit in 
the terminology of the network designer, and 

(b) Numerical integration - the obtaining of a numerically 
accurate solution of tlie initial value probl^ for this 
set of circuit* equations. 

The present state of the art allov/s a reasonably effective solu- 
-ion of the former problem. However, most first generation circuit 
«. axysis prograros are somewhat restricted in scope, as regards the 
latter, in the instance of stiff netvzorks v/hich are characterized 
by widely separated time constants. In such circumstances the 
numerical, integration techniques iioplemonted in tlie first versions 
of programs like SCEPTRE and ECAP require such prohibitively 
small time steps tliat they become impractical as aids to analysis. 

The primary raison d'etre for STICAP is tlie motivation to 
combine tlie capabilities of circuit translation using the state 
variable topological approach v/ith efficient numerical integra- 
tion techniques for transient analysis, using algorithms which 
possess botii stiff and non-stiff capabilities. The STICAP prograra 
is restricted to the analysis of linear time invariant netv/ork.s. 

It represents a merging, with some modifications to each, of 
Pottle's circuit analysis program CORilAP v/ith Gear's program 
ALGO RITE 407 - DIFSUB, for the automatic integration of ordinary 



differential equations . 

The program package is boot viev/ed as consisting of three 
separate components, or nodes of operation, each having some 
advantages and disadvantages over the others in different circum- 
stances. In each mode the common method of circuit translation 
is tlaat originally employed in program CORIIAP; a topological 
approach the result of which is a set of first order linear diffe- 
rential cejuations governing the time evolution of the circuit state 
variajbles. The CORIIAP mode malces selectable the program CORI^AP 
with all previous capabilities, but optional selection of certain 
data printing features. These capabilities include calculation of 
transfer functions, zeroes of tramsmission, frequency and time 
response of tlie circuit. 

The fourth order numerical integration algorithm iirrolemented 
in CORlJhP for time domain analysis is absolutely stable; hence it 
may be used for either stiff or non-stiff iietvrorks. Koi/ever, tiie 
Gtepsize is fixed threugisf nt the duration of computation, a 
feature v;hich can be uneconomical in seme instances. Further, 
other thcin impulse or step functions, the only type of circuit 
input is sampled data. 

The Gear and :.'atrix modes may bo used to compute time domain 
transient, impulse, or step responses only, v;ith tlie option of 
calling CORITAP subroutines to obtain transfer functions and zeroes 
of transmission. The Gear mode allows the selection of either an 
Adam’s integration method, suitable for non- stiff equations, or 
else the methods of Gear, suitable for stiff equations. This mode 
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can be used for analysis of the general linear time invariant net- 
work, v/itli forcing functions specified using the full power of the 
FORTRAN language, or by means of sai'iplcd data. In both cases 
automatic order selection techniques and variations in tiie step 
size are employed as tlie integration precedes, to achieve a desired 
level of accuracy v/ith the minimum niimber of integration steps. 

The maximal order truncation error selectable by changing from one 
algorithm to another via the automatic order selection process is 
an eighth order Adam's method or a sixth order stiff algorithm. 

The matrix method v;as developed by personnel of Old Dominion 
University, Norfolk, Virginia. In tiie matrix mode a spectral decom- 
position of the system matrix in terms of its eigen- values is 
employed, to obt«in a closed form solution which avoids a numerical 
integration. The resulting method is computationally rapid and 
may be used for either stiff or non-stiff netvjorks; hov/ever, it is 
applicedsle only in t!ie case of no repeated eigenvalues of tiie 
system matrix, and for systems v/hose forcing functions are linear 
combinations of sinusoidal, cosinusoids, impulse, or step functions. 

The literature concerning the theoretical aspects of tiie 
combined program package is ratlier extensive; the papers considered 
most helpful are indicated in the bibliographies. In the present 
document tho theori^ first presented by Pottle, Gear, et.al., is to 
some extent paraphrased and/or summarized; in some instances clari- 
fying examples are given. 
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CIRCUIT THEORY 

TliG underlying networl: tlieory and raatliematical algorithitis das- 
cribed by C. Pottlo [10] [11] and incorporated in the design of 
program CORliAP v;ill nou be reiterated, essentially as in Pottle's 
description. The programs used by CORI^AP to tr£mslate the circuit 
description to a set of first order differential equations in the 
state variables of the circuit are coimonly used by all program 
modes. 

2.1 CIRCUIT TRANSLATION: PROM ilSTIJOZ^K DESCRIPTION TO STATE EQUATIONS 
The algorithm for obtaining tlie state equations of a general 
active linear netv/ork innlenented in CORIU\P is a modification of 
tiiat proposed by Dervisoglu [2] and has been described in detail 
elsv/liere [1], [3], [4]. In this iiaplementation it is assumed tliat 
the state equations sought have a sul^set of the capacitor voltages 
and inductor currents as state variables. The problem is, of 
course, to proceed from an easily given description of the network 
to a set of first order coupled differential equations which con- 
tain only state varied:les and independent source inputs. All 
other breincli voltages and currents are removed along the »Jay, ?. 
summary of the metliod follov;s. 

1) Form the fundamental loop matrix with respect to a norraal 
tree; that is, a tree which contains all voltage sources as 
tree branches, all current sources as linl:s, and as many 
capacitive tree branches and inductive linJes as possible. 

This matrix expresses a relation giving all tree-branch 
currents in terms cf link currents and, using its negative 


treuispose> giving all link voltages in terms of tree**branch 
voltages. 

2) Form and solve a set of algebraic equations relating 

a) resistive tree-brainch voltages to their currents (and 
hence to link currents) , 

b) resistive link currents to their voltages (and hence to 
tree-breuich voltages) , 

c) controlled sources to their controlling quantities 
(and hence to link currents and tree-branch voltages) . 

3) Form a set of initial state equations obtained by replacing 
inductor voltages by expressions involving derivatives of 
inductor currents and replacing capacitive currents by expres- 
sions involving derivatives of capacitor voltages. 

4) By row reduction techniques obtain the final state equa- 
tions with dependent state variaibles removed. 

The Fundamental Loop Matrix 

Although a human may find a normal tree by inspection, for a 
computer this tree is roost easily found from the incidence matrix. 
This matrix is easily constructed from the given topological input 
data, each row expressing tlie Current Law with respect to a region 
surrounding one node. A systematic, computer-oriented method for 
obtaining a normal tree and its associated fundamental loop matrix 
from the incidence matrix is as follov;s: 

1) Arrange the columns (branches) of the incidence matrix Q 
in the order 


voltage sources 


capacitors 
resistors 
inductors 
current sources 

Note that since source conversion and transportation are not 
particularly adapted to computerization, the approach taken 
here requires that each branch consist of only orri element; 
voltage and current sources are branches and are ret required 
to be members of "generalized branches" with passive elements. 

2) Apply the standard Gauss- Jordam rov; reduction technique to 
the rows of C v/ith the modification that if, at the i^ step, 
all entries in the i^ column bcloif and including the i^^ rov; 
arc zero, one simply goes on to the first column v/here a 'one' 
can be brought into the i^^ rov7. t^en processing is complete, 
the C matrix is in "row echelon form". 

3) Rearrange the columns (branches) of C to produce a matrix 
of the form 


The first columns witli the single one in them refer to tree 
branches and the rest refer to links. 

Tree branch voltages and link currents will henceforth appear 
in upper case and tree-branch currents and link voltages in lower 
case. With the appearance of the fundamental loop matrix F we have 
obtained the relations 


^ » P ^ and V « -P^ V 


( 1 ) 


relating all lower case quantities to upper case quantities. The 
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extended Gauss- Jordan algorithm appears at other places in the 
sequel and is basic to the v;hole technique. It is introduced here 
without discussion of roundoff errors and such, since at no step 
can numbers other than +1, -1, or 0 appear. The restriction on 
the generality of netvrork v;hich can be handled appearing here is 
that controlled voltage sources must not appear in loops v 'h otiiar 
voltage sources, and controlled current sources must not form cut- 
sets with other current sources. Obviously there are legitimate 
networks excluded by this restriction, but they are rather more 
hypothetical than physical. 

Removal of Resistors and Controlled Sources 

Topological considerations now permit the replacement of any 
lower case quantity by a linear combination of upper case quanti- 
ties. But since each resistor voltage (current) is related to its 
current (voltage) , and each controlled source to its control, a 
set of algebraic equations may be constructed in the form 

- controlled voltage sources 

- tree branch resistors 

- link resistors 

- controlled current sources 

- independent voltage sources 

- tree branch capacitors 

- tree branch inductors (2) 

- link capacitors 

- link inductors 

- independent current sources 
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Figure X Illustrates a case where the first controlled voltage 
source is controlled by the first tree branch capacitori and the 
first link resistor depends (through on (possibly all) the 
tree branch voltages* Only the shaded areas on a given row represent 
possibly non-zero areeus. 


tree 

branches links 



tree 

branches links 



Pig. 1 Formation of Algebraic Eq*-ations 


Once this set of equations hais been formed, it is solved 
using the reduction to roitr echelon fcrti. mentioned above. This 
might also be accomplished using the LU algorithm (5] , which is 
both accurate and efficient. Although in the passive case the 
matrix is positive definite, in the general active case singu- 
larity is a possibility, indicating the fact that one of the quan- 
tities being solved for is arbitrary. Such cases are not likely 
for physical circuits. 

Preliminary State Equations 

The constraint equations for energy storage elements (including 
mutual coupling) are now Introduced: 

* 



( 3 ) 
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The procedure here le to remove ^ and vj, using F and while 
removing dependence on resistive and controlled source quantities 
using the solution of the algebraic equations (2) • and V^. may 
then be removed in favor of state variable derivatives «;)ing the 
matrices €3 and The final results appear in the 



or Sx»/x + Bu 

It should be clear that branch quantities defined its outputs can be 
handled in an identical manner during these reductions, yielding a 
preliminary input-output-state equation 


y-Rx+Cx+Du 

A iw — • a, 


( 5 ) 


F inal State Equations 

The preliminary state equations t^uld be in final form if S 
were lui identity matrix* Since this fact is generally not so, the 
procedure continues by reducing the partitioned matrix fs ! A ! ^ 

Lr* • I iJ 

to row echelon form. Several possibilities nov; arise. 

a) S becomes an identity matrix. In this case S is of full reuik 
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and the final state equations are obtained* 

b) S and A contain at least one sero row* In this case the equation 

of the last row 

b*u-0 (6) 

“n ~ 

implies a dependence among the independent sources. No solu- 
tion exists for such a network* 

c) S but not A contains at least one zero row. In this case the 

last row expresses the equation 

«> “ «n* * + “ (7) . 

is^lying a dependence among the state variables. One may be 
chosen for elimination and the relation above used to remove 
a column from the A and C matrices. Differentiating the 
expression gives 

0 “ Sn^ X.+ E <*) 

Which is used to remove the corresponding element from x and 
hence a column from R and S. Derivatives of input signals may 
thus arise during state variable elimination. XTaen these rela- 
tions have been exiiausted, a smaller matrix IS ( A i Bl is 

Lr I I "j 

again reduced to row echelon form and the process repeated until 
an identity matrix is obtained for S, or all state varleUitles 
disappear. 

In tlieory this reduction technique is a completely general one 
in that the algoritl4A starts with the maximum possible number of 

* 


state variable candidates and can remove as many as required. In 
practice, numerical difficulties can arise, but are entirely con- 
centrated in the row reduction of matrices to rotf echelon form. 

The procedure requires an algorithm for deciding the rank of a 
matrix in the presence of round-off noise. This problem is a funda- 
mental one in numerical analysis which has yet to be solved to the 
satisfaction of all. All network analysis techniques have this 
requirement built in in some form, since in cuiy case the order of 

r 

co^lexity of the notv;ork must be determined. The algorithm under 
discussion here presents the problem in precisely the form which is 
being studied - as a rank determination problem. 

2.2 FROM STATE EQUATIONS TO CRITICAL FREQUENCIES 

The designer of a netxvork euialysis scheme is at this point 
faced with the problem of deciding for just v/hat purposes the pro- 
gram will- be used. Time and frequency response data seem to be the 
most often requested; there are, hov/evcr, users who for one reason 
or another require the critical frequencies of the network. For 
this reason, and because it turns out that it can be done accurately 
without great expenditure of effort, the program proceeds to find 
the natural modes of the network together v;ith the zeros of trans- 
mission between all possible' inputs and outputs. 

One way of locating tlie natural modes of the network is to 
apply an eigenvalue-finding algorithm directly to the A matrix. 
Numerical cmalysts seem to be in rare unanimity that Francis’ QR 
algorithm is the best available [6] , [7] and this algorithm is 
employed by CORNAP. The remaining problem is then to determine the 


zeros. In order to use the QR algorithm for this purpose also, the 
problem of hov7 to construct a matrix whose eigenvalues are tiie zeros 
of transmission of a network must be solved. The cuiswer in this 
c^se lies in considering the "inverse system", a concept which has 
been explored in depth by Brockett [8] • 

The Inverse System 

Consider a network v;ith one input and one output connected in 
the feedback of an operational amplifier of gain k > 0 as shown in 
Fig. 2. 



Pig. 2 Het\ 7 ork in Feedback Loop of Operational Amplifier 
If the system function of tlie original neti^ork is H(s) , the system 
function H(s) of tlie modified network ("inverse system") will be 

H(s) = !S (9) 

1 + kE(s) 

In the limit as k ®, H(s) has the zeros of H(s) as its poles, and 
the poles of H(s) as its zeros. If the A matrix of the inverse sys- 
tem can be found without excessive effprt a solution to our problem 
will have been obtained, since its eigenvalues will be the zeros of 
H(s) we are seeking. Fortunately, the machinery for obtaining this 
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matrix is already available. 


State Equations of the Inverse System 

As can be seen from Fig. 2, from the original preliminary state 
and output-state equations of a single-input/ single-output system 

Sx*=Ax + bu (10) 

T . T 

y*rx + cx + du (11) 


we may by substitution obtain the preliminary state equations of 
the inverse system. Since 


or 


$ = u - k(ft - y) 

A T . T 

= ku - kr X “ k£ “ kdu 

k ,A T. T . 

“ FTkd - £ x) 


( 12 ) 


(13) 


v/e obtain by substitution for u in (10) 


k T • k T V 

S + T r v.^ br X * A - Y .-i ~ v be x + v b u 

^ 1 + kd — — 1 + ]cd — 1 + kd ~ 


(14) 


Before considering how to proceed to the limit k suppose first 

I “ 


that the matrix 


i; ^ 


has by rov; operations been reduced so 


I I . 

that b is a column vector with only one non-zero entry, the last 
(bj^) . The preliminary state equations for the inverse system are 
now no different than for the original, with the exception of the 
last row which becomes 

kb 


m It b 

(s„ 


r^) X 


— n 


n 


k b 

cT) X + ” 


1 + kd - - 1 + kd 


A 

u 


(15) 
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Two cases may now be recognized: 

Cases 1 (d 0) : In passing to the limit k •>, the new last 

row of the state matrices become 



(16) 


Case 2 (d « 0) : Dividing the last roi^ by k and passing to the 

limit# v/e obtain the last rov/ 


S: b^r 

-* n— 

A: -b„c^ (17) 


b: 


b 


n 


Having found preliminary sta+-e equations for the inverse system, we 
may proceed to use the mechanism already at hand f' r producing the 
final state equations with the correct number of ct;ate variables. 

The QR algorithm then finds the zeros of transmission of this parti- 
cular input-output pair. The process must of course be repeated 
for each input-output pair. 


The Gain Constant 

In addition to the poles and zeros of transfer functions which 
have been requested, some sort of gain constant set is required. 
These gain constemts are most e£isily found simultaneously by 


» 
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evaluating 


H(s) 

m0 


C Si - A 


-1 


3 + n 


(18) 


at a single value of s. This value of s has bec^n chosen as real 
(to avoid complex matrix inversion) and positive (to avoid network 
poles) . A search is also made to insure that the value of s chosen 
is not too near a zero of transmission of one of the components of 
H(s). 


Remarks 

In the usual case of proper systems %;ith at most tlie same 
niimber of zeros as poles, the above operations may be performed on 
the final state equations (P-I) with some saving in processing time. 
An improper system with more zeros than poles would in this case 
require a way for the A-matrix to grow in size - no such algorithm 
is included. 

The sequence of events in finding the zeros of a polynomial 
filter (there are, of course, none) is that a zero row appesors at 
the bottom of S at each reduction to row echelon form until at last 
no state variables are left. 

Only the observable, controllable natural modes of a network 
turn as poles of tlie transfer function H(s) . Since the present 
method finds all the natural modes of the netv/ork, we must be 
certain that the uncontrollable or onobservable modes are also 
natural modes of the inverse system, and therefore appear as zeros 
of H(s) to accompany the poles v;hich should not be present. That 
these iiK)des are not altered in going to the inverse system may be 
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seen by considering that the inverse system is obtained by operating 
on the input and output of the original system. Since these modes 
are either uncontrollable or unobservable, they cannot be altered 
by operations on botli input and output. 

2.3 PPEQUENCY M4D TRMJSIENT RESPOHSE CALCUIATIOHS 

Once tlie decision has been made to obtain critical frequencies, 
conventional techniques such as direct evaluation of the transfer 
function at desired frequencies, and partial fraction expemsion 
folloi^ed by evaluation in the time domain, are immediately available 
emd are fairly general. No more general or accurate approach exists, 
in fact, for finding the frequency response, and tlierefore this 
simple approach is the one used in the overall program. A some- 
what more general technique vjas incorporated for computing the 
time response, which is explained in detail elsewhere [9]. The 
method consists of breaking the system function into a cascade of 
second order systems as sliovm in Fig. 3. The state equations of 
each subsystem are solved exactly, vrhile the 



Fig. 3 Transfer Function as a Cascade of Subsyst^=ims 

convolution integral betv/een each subsystem is approximated by a 
fourth-order integration procedure. Sampled input signals are 
easily processed, and tlte method is completely insensitive to 
multiple or clustered poles. 
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CHM>TER III 

STIFFLY STJ^LE IIUIERICAL INTEGRATION 
3.1 INTRODUCTION 

The output from the CORNAP circuit translation routines con- 
sists of tlie differential and algebraic equations of the circuit; 
namely, a state equation 

^ + B^u + Ei§a (3.1) 

and an output equation^ 

y = Cj^x -f Dj^u . (3.2) 

Here 

X is the s X 1 state vector of the circuit? 
u is a T X 1 time varying vector whose components are tlie 
independent sources of the circuit; 
y is a v X 1 vector whose components are the outputs requested 
by the user; and are constant matrices of 

dimensions consistent with the above equations (3.1 - 3.2). 

If is nonzero, the change of variables 

X * q + Ej^u 

transforms the above equations to a form free of source derivatives* 


^It is clear tliat the state cuid output equations for tlie general 
circuit could involve higher derivatives of the inputs; i.e. 

dx — — du d^u 

_ du d^U 

y « C^X + + . . . 

Hov;ever, only equations of the forr.i of (3.1 - 3.2) can be processed 
by STICAP, since no provisions are made to allow user input of 
derivatives of the independent sources. 
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/v3 + Eu 


(3.3) 


Here 


y » Cq + Du • 


(3.4) 


A » C “ 

B * A^Ej^ + D « C^Ej^ + D^, (3.5) 

and q is the transformed state vector. 

The problem of obtaining time domain circuit response is no^.f 
reduced to finding a solution of eons. (3.3 - 3.4). In the Gear 
mode tills is accoitqilished by numerically integrating the initial 
value problem 

^ a Aq + Bu, q(tQ) * Qq* (3.6) 


Stiffly stable numerical integration techniques nay be applied, if 
the system is stiff; othen7ise, an Adam's integration algorithm may 
be selected. In this chapter the necessity of using stiff methods 
is indicated, together v/ith a survey of the related theory. 


3.2 ACCURACY AND STABILITY OF NUIIERICAL HITEGRATIOH 
Consider tlie scalar initial value problem 

= f(t,z), z(tQ) * Zq (3.7) 

consisting of one first order ordinary differential equation, for 
v;hich a solution is desired on the domain t ^ tQ, satisfying ini- 
tially z(tQ) ® Zq. In the absence of a closed form solution one 
settles for some form of a numerically approximate solution 

(y(tj); t^ 1 ^0' ^ 

v;ith the set of time points tj dense enough to yield information 
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sufficient for the purposes at hand. The tj often are chosen uni- 
formly spaced, V7ith « tj^ + h, inhere h is called the step size s. 

The equations employed to generate this approximate solution raay ’3e 
called a numerical integration algorithm . A customary notation for 
the approximating solution values is y(t„) “ y„, t^j « tQ + nh. 

Suppose tlie solution of (3.7) is expressed as 

tn+h 

*n+l * ^n f (u,z(tQ,2Q,u)ldu, (3.£) 

and. let the integral in (3,3) be approxirated by the trapezoidal 
rule 

a+h , 

g(x)dx^^[g(a) +g(a+h)3. (3.9) 

2 

Then we obtain the recursive nuraerical integration algorithm 

yn+l “ yn 

v/here ** f(tn/yn)r ^uid is a numerical approximation to z^. 

Given values yjj, the next solution point, recursively ob- 
tained, is implicitly specified by (3.10). In contrast, the Euler 
algorithr.i 

yn+l “ yn 

e?g?licitly specifies the next solution point. These algorithms pro- 
vide respectively very simple examples of typical members from the 
classes of implicit and explicit linear multistep methods . Funda- 
mental differences in character of these tv/o algorithms ^.-ill now 
be indicated. The phenomena discussed are also characteristic of 


I 


- 21 - 

mor« complicated linear multistep methods.^ 

The accuracy of a numerical Integration algorithm Is customarily 
measured In terms of Its performance when applied to Initial value 
problems (3.7) having as exact solutions the polynomials {l,t^t^r...« 
t^>. Ifr v/hen applied to problems whose solutions are members of 
the test set {l,t,t^,« .. ,t^>, and In the absence of roundoff error 
(l.e., on an Infinite precision machine);, the numerical and the 
exact solutions are Identical » but where error appeairs In the numer- 
ical solution when the exact solution is t^'^’^r the algorithm Is sadd 
to be of order P . Due to the linearity # If the order Is p« the 
linear multistep algorithm, whether Irpliclt or explicit, yields an 
exact solution to the initial value problems v/hose time solution is 
an arbitrary polynomial of degree p or less, in the absence of com- 
putational errors in the initial conditions- 

The order p of the algorithm dictates the magnitude of the 
stepslze h needed to produce a given degree of accuracy, ^dien the 
algorithm Is applied to a problem whose solution Is not a polynomial 
of maximal degree p. This relation may be mere explicitly stated 
by means of the one step truncation error , the error made In comput- 
ing yn+i when exact solution values of all required previous values 
y^-j/ (j “ 0,1,..., k-1; for a k step algorithm) are kno^^n. For 

^The general linear multistep method of k steps is of the form 

j-o 

with the ai,0£ real constants. T.;e method is implicit if 
otherwise, it is explicit. For implicit methods with nonlinear 
fn*f(tnfyn^» some type of predictor corrector iteration must be 
used to solve for y„, given values of (yn-lrYn-a' • • • '^n-k' ^n-l'^n-2' 
...,fj^.)(}. The linear one step methods typified by (3.10) an<" (3.11) 
are considered a subclass of the general family. 
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an algorithm whoso order is p the one step truncation error is of 
the form [ 1 ] 

T - ♦ OOjP^*), (3.12> 

Dhare C- 4 .]^ !• a oonataat and 0 (hP^^) indioataa a function which 
approaches sero as h approaches sero in the same manner as does 

hp+2. 

Thus, the accuracy of the algorithm, as wall as its ecos^resy 
in terms of stepsize magnitude for a given degree of precision, 
increases directly with the integer p. 

It is clear that the es'plicit Euler method (3.11) is of order 
p « 1, and it is readily verified that the trapezoidal method (3.10) 
is of order p « 2 , 

We now examine another characteristic phenomeiia cf these algo- 
rithms; namely, their numerical stability, or as the statistician 


might say, r^ustness under variations in the step size. We shall 
consider as test system the initial value problem 

■ XZ, Z(tQ) * Zq (3.13) 

where X is allo('.*ed to be a real or oomplex number. ^ 


Here the exact i elution is 


X(t-to) , 


z - ZQe""' (3.14) 

the desired sequence to be produced by numerical integration is the 


set of complex numbers 

nXh , . . . 

*n “ *0® ' (tn*to+nh) . 

However, the Euler algorithm produces the sequence 


y„ - a+hX)» Zg, 


^Curiosity concerning the permissiveness of complex values of X 
will be dispelled in the next section. 
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and the trapazoidal rula the sequence 


ti"-' h^2 X>" »o* 


In the limit aa n • and h *► 0 witlt the product nh « ^xT^q 
held ocnatant, both sequences converge *zo the true solution (3.14) 
Boirever, fox h fixed, as n « the Euler sequence diverges, regard- 
less of the size of h, whenever |l i* hX| > 1, while the trapezoidal 
sequence diverges if and only if P?(hX) > 0. This phenomena is 
referred to as numerical instability y when it occurs the accuracy 
of the computation rapidly degenerates. 

For the general linear multistep method, the degree of the 
difference equation usually exceeds that of the differential equa- 
tion it is Intended to approximate. Thiis, certain parasitic 
tions are present, in addition to the desired solution. v?hen ex- 
cited by roundoff or other computational error, these parasitic 
solutions tend to spoil the computation, in proportion to the degree 
that their rate of growth exceeds that of the desired solution. 

A stability requireiT'^t v/hich prevents the spoiling of the com- 
putation by excitation of parasitic solutions will now be deflne<' ; 

A linear multistep method is absolutely stable if the numerical 
solution approaches zero for h fixed, as n approaches infinity, 
applied to all initial value problems (3.13) whose eigenvalue X has 
negati^’e real part; Re(hX) < 0. Should this phenomena occur only 
when hX is restricted to some proper subset D of the left half p2>ine, 
the algorithm is said to be D-strongly stable , or absolutely stable 
on a restricted domain. 


^The numerical solution converges pointwise to the exact solution, 
as h approaches zero. An algorithm characterized by this property 
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The hX plane regions of svahlllty and instability of each al- 
gorithm (3.10) and (3.11) are indicated in Figure 1* 



The trapezoidal rule is absolutely stable/ v/hereas the region of 
stability for the Euler method is the interior of the circle centered 
at C* 0) with radius ^ Euler method is strongly stable 

on the interior of the ci7.cle* For this method the stability re- 
quirements may be expressed in terms of stepsize by the relations 

h < , Re(hX) <■ 0. 

Hence the stepsize h is very restricted by the modulus of the eigen- 
value X of the system being integrated. 

We now make a fundamental observation, exemplified by the pre- 
ceding examples. Namely, for all convergent explicit linear multi- 
step methods, the hX plane stability region is bounded and has the 
origin as either an interior point or a boundary point. This obser- 
vation also holds true for some implicit algorithms, the exception 
being the stiffly stable algorithms of Gear [2] , (31 ; and some 
instants of algorithms of other types, such as the absolutely stable 
Bosenbrock rule reported by Calahan 14]. Hence, for many commonly 
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used techniques (Rungo-Kutta, AcTam's-Bashforth, etc.), there is a 
restriction on the stepsize , to avoid instability, umally of the 
form 

A 

0 < h < j . , A a positive constant , 

where is the largest eigenvalue of the Jacobian matrix of the 
system. In the case of stiff circuits this restriction is peurticu- 
leurly hard to abide, since the step size is determined by the system 
oom^nent with minimum time constant, vrhich in the large contributes 
least to the solution. 

For the general sniff algorithm the stability region is un- 
boimded in the left half plane and possesses a corridor of stable 
approach to the origin (see Pig. 2) . 



Figure 2. Stiff Stability Requireirtents 


In the part of the hX plane to left of the line x ® D, corresponding 
to occurrence of system eigenvalues related to the high frequency 
components, the stiff algoritlim is required to be stable. In the 
region bounded by the lines x =» D, x =* 0 and y » +8, the algorithm 
must be both staJDle and accurate; the rest of the plane is a "Don’t 

* « 
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Care" region and can be stable or unstable. Hence stiff stability 
requires em algorithm which is numerically stcJDle in the region of 
hX corresponding to rapidly decaying system components of diminish- 
ing signif iccince , emd whici: is both accurate and stable in the re- 
gion corresponding to reasonably large step sizes and less slot'/ly 
decaying system components. 

Thus, for physically stable systems, the stability of the stiff 
algorithm is unaffected by the presence of left half plane eigen- 
values of large modulus? for eigenvalues of other types, stability 
can be achieved by adjustments in step size (see Pig. 3 for further 
clarification of this statement) . The off shoot of such methods is 
that the character of the computation is no longer so utterly re- 
stricted by the largest eigenmode; when high frequency effects are 
of diminished importance, the integration can proceed using a time 
step determined only by the needs of accuracy. Economical computa- 
tion times may be achieved by suitably varying the stepsize h and 
the order p (by changing fror. one stiff algorithm to another) as 
the integration proceeds, so as to maximize the step size \ 7 hile 
simultaneously maintaining stability and a desired level of accuracy 
specified by the user [9]. 

The tradeoffs which can be made betvreen the competing factors 
of stability and accuracy, for the stiff algorithms programmed in 
Algorithm 407, are indicated by Figure 3. Typical hX plane stabil- 
ity regions are indicated, for Gear’s algorithms of orders p=2,3,4, 
5,6. In these figures the regions of stability and instability are 
symmetric with respect to the axis of the reals ? only the left half 
plane portions are of any si^j -ficance, assuming physically stable 





Figure 3. Sections of stable regions. 

Curves syiranetric with respect 
to the X-axis. 
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networks. All regions of insted>ility for the lower order xaethods 
are contained within the region of instability of the order p » 6 
algorithm (see Figure 4) • These plots also appear in references 
13), m. 

3.3 TBE NUMERICAL TREATMENT OF A VECTOR SYSTEM 

A study of the previous section presumably arouses curiosity 
concerning the permissiveness of complex X in the test equations 
(3.13) . The dispellation of this curiosity is readily a^ieved when 
one considers the cq>plication of an algorithm for solution of a 
scalar equation to a vector system of simultaneous equations. Ne 
noi-7 provide an exai!^le, again using the trapezoidal rule. 

Suppose the numeric 2 d integration of equation (3.4) 

^ « Aq + Bu, q(tQ) = qQ, 

were to be accomplished by the vector trapezoidal rule, a parallel 
^plication of equation (3.10); namely, 

Jn+l ” ?n + 2<^n+l + » ?o " ?0' (3.15) 

In this instance the s x 1 vectors are defined by 

e AOj^ ••• (3.16) 

This further reduces (3.15) to the form 

(I - jA)Yn+i » (I + 5^)?^ + + %)/ (3.17) 

where I is the identity matrix, s x s. 

Consider for purposes of illustration the case in which the 
eigenvalues X 2 ,X 2 /*«*rXg of A are distinct. Then there exists a 
nonsingular matrix Q for v;hich 
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D a Q“^ AQ a diag(Xi,X 2 r... #Xg) 

Is a diagonal matirix^ v/hose diagonal elements are the eigenvalues 
of A. Then under the transformation 

?n “ 

the algorithm (3.17) becomes 

(1 - 5l»Z„+l - (I + “)2„ + (3.18) 

Thus the transformation uncouples the algorithm, so that there nov? 
arises a system of s simultaneous equations, i « l,2,...,s, 

(1 - |li)Pi,„+i - (1 + (3.19) 

where 

9n ” ? ^i j ,n+”' + • 

Here ’ 

Xj^ is an eigenvalue of A, 

Pi,n is the ith component of z“j^, 

Uj^jj is the jth independent source, evaluated at time n, 

the aj^j are the elements of Q“^B, 2 ind 

the summation on j is over all independent sources. 

Now if (3.19) is to be stable for all bounded independent 
source inputs, it must in particular be stable when these sources 
are switched off; i.e., the natural modes of the algorithm must be 
stable. Hov/ever, setting gjj * 0 in (3.19), it emerges that a 
necessity 2 uid sufficient condition for stability is that for each 
i “ 1,2, 3,..., 5, Re(hXj^) be negative. But this is the same criterion 
as was estciblished for the scalar case. Moreover, the result is 
equally valid for occurrence of repeated eigenvalues, as can readily 
be established. 
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3.4 THE STIFF /ALGORITHMS OF GE/AR 

The general form of the vector implicit linear multistep algor- 
ithm for solution of equation (3.6), 

^ « Aq + Bu, q(to) » qQ, (3.20) 

is given by tlie relation ^ 

®K^n+k ®k-iyn+k-l “ ^t^k^n+k • •'•’^o^nJ • (3.21) 

Here a]^3]c ^ Or the functions fn-t>j satisfy^ 

^n+j “ ^^n+j j ' 3 ** 0/1^2, ...,k 

and yji is the numerical approximation of q^* 

The algorithm is of order p if the relations Cq^Cj^*. , ,ssCp=0, 

^ 0, are satisfied, v/here [1] 

Cq * Oq ■** ®1 '*' • • • ®k 

^1 “ ®1 *** ^^2 + • • • + + Ox + ••• + Ok) 

(3.22) 

Cq ** ^(ttx+2^a2+. . .+k*3a3^) - (q-i) |' f0x+2‘J“^S2+« • .+k^“^3j^l . 

Once having restricted the a*s and 3's to produce an algorithm 
of order p, it emerges that the B values are uniquely specified in 
terms of the a values by eqns. (3.22) . Conversely r the additional 
specification tliat p *= k allows a unique specification of the a 
values in terms of the B values. In either case, for a fixed order 


^In this section the theory of linear multistep algorithms, special 
izcd to the case of a linear dynamical system, is given. The gen- 
eral theory, as indicated in the literature tabulated in the bibli 
ography, applies equally v;ell to the nonlinear case. 
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p, the arbitrary parameters may then be chosen to obtain satisfactory 
stability properties. Criteria upon which such a value judgment may 
be based will* noiir be indicated. 

The numerical stability of tlie algorithm, as defined in the 

previous sections, is governed by the follov/ing theorem [2] [5]: 

\ * 

Stability Theorem 

The numerical solution of (3.21) is stable if and only if tlia 
stepsize h is restricted so that the roots of the polynomial equa- 
tions 

k . 

Z - hXiBk_^y''■3 - 0, (3.23) 

j=0 ■’ i 

(v;ith Ai, i = 1,2,... ,s, the eigenvalues of the matrix h) are within 
the unit circle in the complex plane, or else on the unit circle and 
not repeated. 

In view of the foregoing discussion, it appears tliat in order 
to design the optimal k-step algorithm one should maximize simul- 
taneously the integer p (for accuracy) and the region of the complex 
hX pleuie which yields a stable root configuration for each of eqns. 
(3.23). Ho’./ever, a theorem of Dahlquist (1] concerned with con- 
vergence restricts the maximal p, for k fixed, to either k + 1 or 
k + 2, depending upon whether k is cdd or even. Fixing p at its 
maximal value, one next seeks the parameter choice yielding the 
broadest range of stability. Assuming a physically stable differen- 
tial equation, absolute stability of the difference equation is the 
natural criterion, for broad general application. Unfortunately, 
again a result of Dalilquist [5] limits k to the value k « 2, p £ 2, 
for absolute stability. Hov/ever, the loi7 order involved here implies 
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an algorithm of limited usefulness. Thus tlie motivation for dis- 
covery of the stiff metliods arises; methods not quite absolutely 
stable, but strongly stable on a restricted domain of the left half 
plane, xvith the characteristic corridor of stable approach to the 
origin (Fig. 2). 

For stiff algorithms, a result of TJidland [6] imposes the re- 
striction p £ k, for k > 1; p * 2, for I: = 1. Specific stiff alg.v:: 
ithms with p = k, k = 1 (1) 8 have been discovered oy Gear [3], (7j ; 
the general question is still tJie subject of research. 

The algorithms discovered by Gear, v;ith p » k; k = 1,2, 3, 4, 5, 6 
are implemented in ALGORITHIi 407. These algorithms are of the form 

^n+l * ^‘k-l yn + “k-2 yn-1 +•••+ «oyntl-k ^®k^n+l' (3.24) 

and require, in the programmed formulation (7) , the fevjest function 
evaluations and saving of least inforr^ation from step to step of all 
stiff algorithms in the class p * k. The values of the a's and B's 
are given in Table I. A sample plot of the stability region, for 
k » 6, is given in Fig. (4). The stability regions for each algor- 
ithm appear in Fig. (3) (the stability plots are symmetric vjith 
respect to tie real axis). 

Table I. Coefficients of Stiff Algorithms 


K 

i ^0 

“1 j 

»2 

“3 

«4 

“5 1 


2 

2/3 

4/3 

-1/3 

0 

0 

0 

0 

1 

3 

6/11 

10/11 

-9/11 

2/11 

0 

1 

0 

0 

4 

12/25 

48/25 

-36/25 

16/25 

-3/25 

0 

0 

5 

60/137 

300/137 

-300/137 

200/137 

-75/137 

12/137 

0 

6 

60/147 

360/147 , 

-450/147 

400/147 

-225/147 

72/147 

-10/147 
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These algorithms are obtained as follows: Equations (3,23) 



ajcP*^ + + ttg 

+ • • • + $Q 


(3.25) 


of the complex root plane (the p plane) onto the H * hX plane. i\s 
h approaches infinity the roots of (3.23) approach the poles of H; 


as h approaches zero the roots of (3.23) approach the zeroes of I!. 
Since the roots of a polynomial vary continuously v;ith its coeffi- 
cients, the choice of tlie poles of H to be strictly within the unit 


circle assures that each hX point in some neighborhood of infinity 


will be stable. By continvdty, this stability region can be extended 
so that at least each point exterior to the image of the unit circle 
|p| » 1 is a stedsle point. (Moreover, v/hen H is one-to-one on the 
unit circle v?ith poles inside this image locus separates the stable 
cuid unstable regions [83 of the hX plane.) The unique stability 
mapping obtained by choosing all poles of II at the origin and deter- 
mining the zeroes by the req* Bremen t that the order of the corre- 
sponding algorithm be p » k yields stiffly stable algorithi.is , for 
k = 1,2, 3, 4, 5, 6, characterized by the stability regions of Figure 4. 


3.5 PROGRAMIIED FORT lULAT ION 

In the preceding sections v/e have presented the general theory 
underlying the linear multistep integration techniques exploited by 
Gear, and indicated the basic algorithms implemented in ALGORITHM 
407. The mathematical formulation of these stiff methods and the 
Adam's- Bash forth methods for non-stiff equations actually programmed 
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will not be discussed here>^ The Interested reader may refer to 
references [9], [10] for a detailed discussion of the mathematical 
formulation employed, the error controls used, and criterion for 
varying the stepsise and order of the metliods. 

As indicated in the preceding discussion, the stiff system 
theory applies equally vfell to nonlinear systems, and ALGORITHM 407 
is programmed for sudi. Hence it might be argued that the predictor * 
corrector techniques used for solving equation (3.21) iteratively 
for yjj+ 3 ^ in terms of preceding values could perhaps be inefficient 
when applied to a linear system. However this conjecture is not 
justified; in the linear case tiae Newton iteration involved converges 
with only one iteration [9]. Hence the process reduces to a matrix 
inversion at each step, which cannot be avoided, regeurdless of the 
linearity, if a linear multi step method is assumed. Furthermore, 
in its programmed version the mathematical formulation automatically 
provides the infonaation n<'.eded for determining the necessity of a 
change in stepsize or order of the method, an essential feature of 
€iny numerical integration technique which is to be effective over 
a broad range of netifjork analysis problems. 


%e remark that the formulation employed is chosen for facility 
with vjhicli it lends itself to stepsize and order changing. 
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CHAPTER IV 
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THE ilATRIX METHOD 

4.1 INTRODUCTION 

Aa indicated in Chapter 3, the differential and algebraic equa- 
tions of a network processable by STICAP may be described in tlie 
form 

dx 

- Ax t Bu, x(tQ> • Xq (4.1) 

y ■ cx - Du, 


where x is the state vector (or transformed state vector) ; y is the 
vector of outputs; u is tlie vector of independent sources; A, B, C, 
D are oonstemt matrices. 

The solution of the state equation is 


t 

X • (exp(A(t-to) } Xq + / exp(A(t - t) v(t) dx, (4.2) 

to 

where 

v(t) » Bu(t) . 

One way to avoid the problems v/hich arise in numerical integration 
of stiff circuits is to integrate (4.2) exactly for certain wave- 
forms such as sinusoids and block pulses, which constitute many of 
the excitations encountered in electronic circuits. 

Tills explicit integration ceui be carried out if the transition 
matrix exp(A(t - t) Ccui be expressed as a linear combination of the 
eigenmodes of the system. The eigenvalues of the A matrix are 
obtained by tlie QR transformation, the best method available. The 
particular expansion in tcn«is of the eigenmodes essential for the 
STICAP matrix algorithm was studied by Kirchner [1]. In tlie present 
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algorithm the axpanaion is limited to the case of non-degonerate 
modes; the general case is more complexi but still solvable. 

4.2 SXG&mtODE BXPANSI(»I OF exp (At) 

Consider tho case for which the e:i jenvalues of A are distinct. 
For this case the transition matrix is ci^ressible as 


n n 
exp(At) • y I 

i*l j«l 


Ci^A 


exp(Xit) 


(4.3) 


The coefficients form a matrix which turns out to be the inverse 
of the Vandermonde matrix 


C » CC^j) 




(4.4) 


Tlie Vandermonde inverse can be confuted efficiently through the 
method developed by Kaufman (2} , 








n-j-k 


k»0 


ir (X^ -Xj^) 
k-1 


k;<l 


(4.5) 


where the aj^ values arc the coefficients of the characteristic 
equation 


P(X) » a^X*' + aiX*^-l + . . . + a^j.j^X + a^j (4.6) 

with Sq » 1 



39 


Tho aj^ coefficients can be formed from the combination of the 
traces of A , k = 1, . . . / n 


^ * m=l ^ ® ® 


( 4 , 7 ) 


- trA® = i (X®) 
£«1 *- 


( 4 . 8 ) 


The total number of operations involved in forming the C matrix is 
about 4.5 operations, which is considerably faster them Gauss 
elimination, for large N. It Is very cumbersome to store the transi- 
tion matrix as expressed by Equation (4.3). Since the matrix is 
always multiplied into a vector such as exp(A(t-T) )y, then it would 
be convenient to consider the vectors 

exp(A(t-T))v * £^exp(X^(t-T)) + . . .+ (4.S) 

where 


Ej “ <Cjl^ + CjjA + • • • + CjnA""^)v 


(CjiYi 


r + . . • + 


'j2*2 


Dn*n' 


(4.10) 


with the Yift+i vector difined by 


Yl = V 


^m+1 " ^ ^m-1 


(4.11) 


This can be easily set up as a recursive process, if use is made to 
take advantage of the sparseness of the A matrix so that the program- 
ming is more efficient. 



The integral in Equation (4.2) can no»7 be expressed as 


t n t 

/ exp{A(t-x)) V (x) dx « r / p. exp(X.(t-x)) dx (4.12) 

to 1-1 to 3 3 

For a finite number of excitation sources pj can be written as 


Ic«M 

Pj “ tkS„(T) (4.13) 

ks»l 

where Fj- is a constant coluion vector, and Equation (4,12) becomes 


n M t 

Z E / expUj (t“X) )sj^(x) dx (4,14) 

j*l k*l tQ 


Notice the integrals in Equation (4.14) are now scalar quantities and 
these are nothing more than convolution integrals. Since the inte- 
gral 

t t to 

/f dx = / f dx - / f dx (4.15) 

to o o 


it is sufficient to examine the integral 


t 

/f expUj (t-x) )sj^(x) dx 
to 


(4.16) 


from which Equation (4,15) can be computed. The Laplace treinstorm of 
the above integral is 


_ L- _ !L- b(P) 

p-Xj ® p-Xj d(p) “ p-Xj ■*' d(p) 


(4.17) 


The first terra rn the above equation is the natural mode in the cir- 
cuit while the remainder term is the response due to tl'^s forcing 
function. If Xj and the poles of the driving function are situated 
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widely apart, er. the ill-conditioned case comes up* In the 
present algor: . an exact explicit integration is performed for 

Equation 9*16)r for certain excitation waveforms, thereby eliminating 
the difficulty. I^hen there are complex eigenvalues, t!ien tlie 
following integral pair is considered 

t ^ * A 

{/ exp(Xi(t-T) ) rs (t) dt + / exp<Xj^(t-T) )r*s(T) dx} (4.1Q,» 

o o 

where r is an element in the column vector of Equation {4*14). 
Equation (4*10) is valid if the original state equations consist 
of real quantities only. For a full discussion of the algorithm 
mathematics concerning the exact integration of (4.14) for the 
functions STICAP allowable as excitations see [31, [4], [5]. 
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